/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
** ****************************************************************** */

//
// 文件: GPUElementKernels.cu
// 创建时间: 2025年10月2日
// 说明: GPU单元刚度矩阵计算核函数实现
// 作者: OpenSees GPU加速项目组
//
// 功能:
//   - 实现GPU单元刚度矩阵计算核函数
//   - 支持Truss、Beam、Quad等常用单元类型
//   - 提供高性能并行刚度矩阵计算
//

// ============================================================
// Performance optimization: Disable kernel debug output
// Include GPU_ENABLE_DEBUG_OUTPUT from GPUSolver.h
// ============================================================
#ifndef GPU_ENABLE_DEBUG_OUTPUT
#define GPU_ENABLE_DEBUG_OUTPUT 0
#endif

#if GPU_ENABLE_DEBUG_OUTPUT
#define KERNEL_DEBUG_PRINT(...) KERNEL_DEBUG_PRINT(__VA_ARGS__)
#else
#define KERNEL_DEBUG_PRINT(...) ((void)0)
#endif

#include "GPUElementKernels.h"
#include "GPUConcreteKernels.h"
#include <cuda_runtime.h>
#include <stdio.h>

//==============================================================================
// GPU辅助设备函数实现
//==============================================================================

/**
 * @brief 计算平面应力材料矩阵
 *
 * 平面应力假定：σz = τxz = τyz = 0
 * 材料矩阵D为3x3，用于平面问题
 */
__device__ void computePlaneStressD(
    double E,      // 弹性模量
    double nu,     // 泊松比
    double D[3][3] // 输出：材料矩阵
)
{
    // 平面应力材料矩阵系数
    double factor = E / (1.0 - nu * nu);

    // 材料矩阵D（对称矩阵）
    D[0][0] = factor;              // E/(1-ν²)
    D[0][1] = factor * nu;         // νE/(1-ν²)
    D[0][2] = 0.0;

    D[1][0] = factor * nu;         // νE/(1-ν²)
    D[1][1] = factor;              // E/(1-ν²)
    D[1][2] = 0.0;

    D[2][0] = 0.0;
    D[2][1] = 0.0;
    D[2][2] = factor * (1.0 - nu) / 2.0;  // G = E/(2(1+ν))
}

/**
 * @brief 从3D切线刚度矩阵(6x6)提取平面应力刚度矩阵(3x3)
 *
 * 使用静力凝聚消除out-of-plane分量 (εzz, γxz, γyz)
 * D_plane = D11 - D13 * inv(D33) * D31
 *
 * @param D_3d 输入：3D切线刚度矩阵 [6x6]，顺序：εxx,εyy,εzz,γxy,γyz,γxz
 * @param D_plane 输出：平面应力刚度矩阵 [3x3]，顺序：εxx,εyy,γxy
 */
__device__ void staticCondensation3Dto2D(
    const double D_3d[36],  // 输入：6x6矩阵（行优先）
    double D_plane[3][3]    // 输出：3x3矩阵
)
{
    // OpenSees应变顺序：εxx(0), εyy(1), εzz(2), γxy(3), γyz(4), γxz(5)
    // 平面应力：保留 εxx, εyy, γxy；消除 εzz, γyz, γxz

    // 提取D11 (3x3)：[εxx,εyy,γxy] 对应索引 [0,1,3]
    double D11[3][3];
    int idx_2d[3] = {0, 1, 3};  // εxx, εyy, γxy在6D中的索引
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            D11[i][j] = D_3d[idx_2d[i] * 6 + idx_2d[j]];
        }
    }

    // 提取D13 (3x3)：[εxx,εyy,γxy] vs [εzz,γyz,γxz]
    double D13[3][3];
    int idx_3d[3] = {2, 4, 5};  // εzz, γyz, γxz在6D中的索引
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            D13[i][j] = D_3d[idx_2d[i] * 6 + idx_3d[j]];
        }
    }

    // 提取D33 (3x3)：[εzz,γyz,γxz] vs [εzz,γyz,γxz]
    double D33[3][3];
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            D33[i][j] = D_3d[idx_3d[i] * 6 + idx_3d[j]];
        }
    }

    // 计算 inv(D33) - 使用3x3矩阵求逆
    double det = D33[0][0] * (D33[1][1]*D33[2][2] - D33[1][2]*D33[2][1])
               - D33[0][1] * (D33[1][0]*D33[2][2] - D33[1][2]*D33[2][0])
               + D33[0][2] * (D33[1][0]*D33[2][1] - D33[1][1]*D33[2][0]);

    if (fabs(det) < 1e-20) {
        // D33奇异，回退到D11（相当于忽略out-of-plane耦合）
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                D_plane[i][j] = D11[i][j];
            }
        }
        return;
    }

    double inv_D33[3][3];
    inv_D33[0][0] = (D33[1][1]*D33[2][2] - D33[1][2]*D33[2][1]) / det;
    inv_D33[0][1] = (D33[0][2]*D33[2][1] - D33[0][1]*D33[2][2]) / det;
    inv_D33[0][2] = (D33[0][1]*D33[1][2] - D33[0][2]*D33[1][1]) / det;
    inv_D33[1][0] = (D33[1][2]*D33[2][0] - D33[1][0]*D33[2][2]) / det;
    inv_D33[1][1] = (D33[0][0]*D33[2][2] - D33[0][2]*D33[2][0]) / det;
    inv_D33[1][2] = (D33[0][2]*D33[1][0] - D33[0][0]*D33[1][2]) / det;
    inv_D33[2][0] = (D33[1][0]*D33[2][1] - D33[1][1]*D33[2][0]) / det;
    inv_D33[2][1] = (D33[0][1]*D33[2][0] - D33[0][0]*D33[2][1]) / det;
    inv_D33[2][2] = (D33[0][0]*D33[1][1] - D33[0][1]*D33[1][0]) / det;

    // 计算 D13 * inv(D33) (3x3矩阵)
    double D13_invD33[3][3];
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            D13_invD33[i][j] = 0.0;
            for (int k = 0; k < 3; k++) {
                D13_invD33[i][j] += D13[i][k] * inv_D33[k][j];
            }
        }
    }

    // 计算 D13 * inv(D33) * D31 (注意D31 = D13^T)
    double correction[3][3];
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            correction[i][j] = 0.0;
            for (int k = 0; k < 3; k++) {
                correction[i][j] += D13_invD33[i][k] * D13[j][k];  // D13[j][k] = D31[k][j]
            }
        }
    }

    // 最终结果：D_plane = D11 - correction
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            D_plane[i][j] = D11[i][j] - correction[i][j];
        }
    }
}

/**
 * @brief 计算2D四边形单元的Jacobian矩阵
 *
 * Jacobian矩阵将自然坐标系(ξ,η)映射到物理坐标系(x,y)
 * J = [∂x/∂ξ  ∂y/∂ξ ]
 *     [∂x/∂η  ∂y/∂η ]
 */
__device__ int computeQuadJacobian(
    double coords[4][2],  // 输入：单元节点坐标
    double xi,            // 输入：自然坐标ξ
    double eta,           // 输入：自然坐标η
    double J[2][2],       // 输出：Jacobian矩阵
    double& detJ          // 输出：Jacobian行列式
)
{
    // 形函数对自然坐标的导数
    double dN_dxi[4], dN_deta[4];

    dN_dxi[0] = -0.25 * (1.0 - eta);
    dN_dxi[1] =  0.25 * (1.0 - eta);
    dN_dxi[2] =  0.25 * (1.0 + eta);
    dN_dxi[3] = -0.25 * (1.0 + eta);

    dN_deta[0] = -0.25 * (1.0 - xi);
    dN_deta[1] = -0.25 * (1.0 + xi);
    dN_deta[2] =  0.25 * (1.0 + xi);
    dN_deta[3] =  0.25 * (1.0 - xi);

    // 计算Jacobian矩阵
    J[0][0] = 0.0;  // ∂x/∂ξ
    J[0][1] = 0.0;  // ∂y/∂ξ
    J[1][0] = 0.0;  // ∂x/∂η
    J[1][1] = 0.0;  // ∂y/∂η

    for (int i = 0; i < 4; i++) {
        J[0][0] += dN_dxi[i]  * coords[i][0];
        J[0][1] += dN_dxi[i]  * coords[i][1];
        J[1][0] += dN_deta[i] * coords[i][0];
        J[1][1] += dN_deta[i] * coords[i][1];
    }

    // 计算Jacobian行列式
    detJ = J[0][0] * J[1][1] - J[0][1] * J[1][0];

    // 检查奇异性
    if (fabs(detJ) < 1e-14) {
        return -1;  // 奇异Jacobian矩阵
    }

    return 0;  // 成功
}

/**
 * @brief 计算2D四边形单元的形函数导数（物理坐标系）
 *
 * 通过Jacobian逆矩阵将自然坐标系导数转换为物理坐标系导数
 * [∂N/∂x] = J⁻¹ [∂N/∂ξ ]
 * [∂N/∂y]        [∂N/∂η ]
 */
__device__ void computeQuadShapeFunctionDerivatives(
    double xi,            // 输入：自然坐标ξ
    double eta,           // 输入：自然坐标η
    double invJ[2][2],    // 输入：Jacobian逆矩阵
    double dN_dx[4],      // 输出：形函数对x的导数
    double dN_dy[4]       // 输出：形函数对y的导数
)
{
    // 形函数对自然坐标的导数
    double dN_dxi[4], dN_deta[4];

    dN_dxi[0] = -0.25 * (1.0 - eta);
    dN_dxi[1] =  0.25 * (1.0 - eta);
    dN_dxi[2] =  0.25 * (1.0 + eta);
    dN_dxi[3] = -0.25 * (1.0 + eta);

    dN_deta[0] = -0.25 * (1.0 - xi);
    dN_deta[1] = -0.25 * (1.0 + xi);
    dN_deta[2] =  0.25 * (1.0 + xi);
    dN_deta[3] =  0.25 * (1.0 - xi);

    // 转换到物理坐标系
    for (int i = 0; i < 4; i++) {
        dN_dx[i] = invJ[0][0] * dN_dxi[i] + invJ[0][1] * dN_deta[i];
        dN_dy[i] = invJ[1][0] * dN_dxi[i] + invJ[1][1] * dN_deta[i];
    }
}

/**
 * @brief 计算应变-位移矩阵（B矩阵）
 *
 * B矩阵将节点位移映射到应变：ε = B * u
 * 对于2D问题（4节点，8DOF）：
 * B = [dN₁/dx    0      dN₂/dx    0      dN₃/dx    0      dN₄/dx    0   ]
 *     [   0    dN₁/dy     0    dN₂/dy     0    dN₃/dy     0    dN₄/dy]
 *     [dN₁/dy  dN₁/dx  dN₂/dy  dN₂/dx  dN₃/dy  dN₃/dx  dN₄/dy  dN₄/dx]
 */
__device__ void computeStrainDisplacementMatrix(
    double dN_dx[4],  // 输入：形函数对x的导数
    double dN_dy[4],  // 输入：形函数对y的导数
    double B[3][8]    // 输出：应变-位移矩阵
)
{
    for (int i = 0; i < 4; i++) {
        int col = 2 * i;

        // εxx = ∂u/∂x
        B[0][col]     = dN_dx[i];
        B[0][col + 1] = 0.0;

        // εyy = ∂v/∂y
        B[1][col]     = 0.0;
        B[1][col + 1] = dN_dy[i];

        // γxy = ∂u/∂y + ∂v/∂x
        B[2][col]     = dN_dy[i];
        B[2][col + 1] = dN_dx[i];
    }
}

/**
 * @brief 从全局位移向量中提取单元节点位移
 *
 * @param d_nodeDisplacements 全局节点位移数组（所有DOF）
 * @param DOFIndices 单元DOF索引数组（单元局部DOF到全局DOF的映射）
 * @param numDOF 单元DOF数量
 * @param elementDisplacements 输出：单元节点位移
 */
__device__ void extractElementDisplacements(
    const double* d_nodeDisplacements,
    const int* DOFIndices,
    int numDOF,
    double* elementDisplacements
)
{
    if (d_nodeDisplacements == nullptr) {
        // 如果位移数据不可用，设置为零
        for (int i = 0; i < numDOF; i++) {
            elementDisplacements[i] = 0.0;
        }
        return;
    }

    // 从全局位移向量提取单元位移
    for (int i = 0; i < numDOF; i++) {
        int globalDOF = DOFIndices[i];
        if (globalDOF >= 0) {
            elementDisplacements[i] = d_nodeDisplacements[globalDOF];
        } else {
            elementDisplacements[i] = 0.0;
        }
    }
}

/**
 * @brief 使用B矩阵计算应变：ε = B * u
 *
 * 对于平面应力问题：
 * - 应变向量：[εxx, εyy, γxy]（3分量）
 * - B矩阵：3x8（对于4节点单元）
 * - 位移向量：[u1, v1, u2, v2, u3, v3, u4, v4]（8分量）
 *
 * 对于3D本构需要扩展为6分量：[εxx, εyy, εzz, γxy, γxz, γyz]
 * 其中 εzz = γxz = γyz = 0（平面应力条件）
 *
 * @param B 应变-位移矩阵（3x8）
 * @param elementDisplacements 单元节点位移（8分量）
 * @param strain 输出：应变向量（6分量，用于3D本构）
 */
__device__ void computeStrainFromBMatrix(
    double B[3][8],
    const double* elementDisplacements,
    double* strain,
    double nu  // 添加泊松比参数用于计算εzz
)
{
    // 计算平面应变分量：ε = B * u
    strain[0] = 0.0;  // εxx
    strain[1] = 0.0;  // εyy
    strain[3] = 0.0;  // γxy

    for (int i = 0; i < 8; i++) {
        strain[0] += B[0][i] * elementDisplacements[i];  // εxx
        strain[1] += B[1][i] * elementDisplacements[i];  // εyy
        strain[3] += B[2][i] * elementDisplacements[i];  // γxy
    }

    // ✅ 平面应力条件修正：求解εzz使得σzz=0
    // 对于弹性材料的平面应力：εzz = -ν/(1-ν) * (εxx + εyy)
    // 这是PlaneStressMaterial包装器通过Newton迭代得到的结果
    strain[2] = -nu / (1.0 - nu) * (strain[0] + strain[1]);  // εzz (平面应力)
    strain[4] = 0.0;  // γxz = 0 (平面应力假设)
    strain[5] = 0.0;  // γyz = 0 (平面应力假设)
}

//==============================================================================
// GPU单元刚度矩阵计算核函数实现
//==============================================================================

/**
 * @brief GPU并行计算四边形单元刚度矩阵
 *
 * 核心算法：
 * 1. 提取单元几何和材料参数
 * 2. 2x2高斯积分点循环
 * 3. 每个积分点：计算Jacobian、B矩阵、K += B^T*D*B*detJ*t*w
 * 4. 将单元刚度矩阵写入全局缓冲区
 *
 * 性能优化：
 * - 使用共享内存存储单元刚度矩阵
 * - 每个线程处理一个单元
 * - 避免分支divergence
 */
__global__ void computeQuadStiffness(
    GPUElementData* d_elementData,              // 输入：单元数据数组
    double* d_stiffness,                        // 输出：单元刚度矩阵缓冲区
    double* d_nodeDisplacements,                // 输入：节点位移（几何非线性，可选）
    ConcreteState* d_materialStates_committed,  // 输入：上次收敛的材料状态（只读）
    ConcreteState* d_materialStates_trial,      // 输出：当前试算的材料状态（写入）
    double* d_strainTotal,                      // 输入：上次收敛的应变（只读）
    double* d_strainCurrent_trial,              // 输出：当前试算应变（写入）
    int numElements,                            // 输入：单元数量
    int maxDOFPerElement                        // 输入：缓冲区布局的最大DOF数
)
{
    // 线程索引 = 单元索引
    int elemIdx = blockIdx.x * blockDim.x + threadIdx.x;
    if (elemIdx >= numElements) return;

    // 获取当前单元数据
    GPUElementData& elem = d_elementData[elemIdx];

    // ====== 调试：检查是否进入kernel ======
    if (elemIdx == 0) {
        KERNEL_DEBUG_PRINT("[DEBUG] computeQuadStiffness: elemIdx=0 entered\n");
        KERNEL_DEBUG_PRINT("  d_elementData pointer: %p\n", d_elementData);
        KERNEL_DEBUG_PRINT("  elem address: %p\n", &elem);
        KERNEL_DEBUG_PRINT("  elem.elementTag=%d\n", elem.elementTag);
        KERNEL_DEBUG_PRINT("  elem.elementType=%d\n", elem.elementType);
        KERNEL_DEBUG_PRINT("  elem.numNodes=%d\n", elem.numNodes);
        KERNEL_DEBUG_PRINT("  elem.numDOF=%d\n", elem.numDOF);
    }

    // 检查单元类型是否为Quad
    GPUSupportedElementType elemType = getGPUElementType(elem.elementType);

    // ====== 调试：输出类型转换结果 ======
    if (elemIdx == 0) {
        KERNEL_DEBUG_PRINT("[DEBUG] computeQuadStiffness: elemIdx=0, GPUSupportedElementType=%d\n", (int)elemType);
        KERNEL_DEBUG_PRINT("[DEBUG] Expected types: FOUR_NODE_QUAD=%d, FOUR_NODE_QUAD_3D=%d, FOUR_NODE_QUAD_02=%d\n",
               (int)GPUSupportedElementType::FOUR_NODE_QUAD,
               (int)GPUSupportedElementType::FOUR_NODE_QUAD_3D,
               (int)GPUSupportedElementType::FOUR_NODE_QUAD_02);
    }

    if (elemType != GPUSupportedElementType::FOUR_NODE_QUAD &&
        elemType != GPUSupportedElementType::FOUR_NODE_QUAD_3D &&
        elemType != GPUSupportedElementType::FOUR_NODE_QUAD_02) {
        if (elemIdx == 0) {
            KERNEL_DEBUG_PRINT("[DEBUG] computeQuadStiffness: elemIdx=0 SKIPPED (type mismatch)\n");
        }
        return;  // 不是Quad单元，跳过
    }

    if (elemIdx == 0) {
        KERNEL_DEBUG_PRINT("[DEBUG] computeQuadStiffness: elemIdx=0 PASSED type check\n");
    }

    // ====== 步骤1：提取节点坐标（4个节点） ======
    double coords[4][2];
    for (int i = 0; i < 4; i++) {
        coords[i][0] = elem.nodeCoords[i * 3 + 0];  // x坐标
        coords[i][1] = elem.nodeCoords[i * 3 + 1];  // y坐标
    }

    // ====== 步骤2：提取材料参数 ======
    double E  = elem.materialProps[0];   // 弹性模量
    double nu = elem.materialProps[1];   // 泊松比
    double t  = elem.materialProps[8];   // ✅ 厚度（修复：从[8]读取，不是[2]）

    // ====== 调试：检查材料参数 ======
    if (elemIdx == 0 && threadIdx.x == 0) {
        KERNEL_DEBUG_PRINT("[GPU_MAT] elem=%d: E=%.6e nu=%.6e t=%.6e\n",
               elemIdx, E, nu, t);
    }

    // ====== 步骤3：计算平面应力材料矩阵 ======
    double D[3][3];
    computePlaneStressD(E, nu, D);

    // ====== 步骤3.5：提取单元节点位移（用于非线性材料） ======
    double elementDisplacements[8];
    extractElementDisplacements(d_nodeDisplacements, elem.DOFIndices, 8, elementDisplacements);

    // 调试输出（仅第一个单元）
    if (elemIdx == 0 && threadIdx.x == 0) {
        KERNEL_DEBUG_PRINT("[GPU_DISPL] elem=%d: DOFIndices=[%d,%d,%d,%d,%d,%d,%d,%d]\n",
               elemIdx, elem.DOFIndices[0], elem.DOFIndices[1], elem.DOFIndices[2], elem.DOFIndices[3],
               elem.DOFIndices[4], elem.DOFIndices[5], elem.DOFIndices[6], elem.DOFIndices[7]);
        KERNEL_DEBUG_PRINT("[GPU_DISPL] elem=%d: u=[%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e]\n",
               elemIdx, elementDisplacements[0], elementDisplacements[1], elementDisplacements[2], elementDisplacements[3],
               elementDisplacements[4], elementDisplacements[5], elementDisplacements[6], elementDisplacements[7]);
    }

    // ====== 步骤4：使用共享内存存储单元刚度矩阵 ======
    // 注意：8x8矩阵需要64个double，对于大block size可能超出共享内存限制
    // 这里使用线程局部内存，性能略低但更安全
    double K_elem[8][8];

    // 初始化单元刚度矩阵为0
    for (int i = 0; i < 8; i++) {
        for (int j = 0; j < 8; j++) {
            K_elem[i][j] = 0.0;
        }
    }

    // ====== 步骤5：2x2高斯积分点循环 ======
    // 高斯积分点和权重
    const double gp = 0.577350269189626;  // 1/√3
    double xi_vals[4]  = {-gp,  gp,  gp, -gp};
    double eta_vals[4] = {-gp, -gp,  gp,  gp};
    double weights[4]  = {1.0, 1.0, 1.0, 1.0};

    // 遍历4个积分点
    for (int gp_idx = 0; gp_idx < 4; gp_idx++) {
        double xi  = xi_vals[gp_idx];
        double eta = eta_vals[gp_idx];
        double w   = weights[gp_idx];

        // ====== 步骤5.1：计算Jacobian矩阵 ======
        double J[2][2], detJ;
        int jacResult = computeQuadJacobian(coords, xi, eta, J, detJ);

        if (jacResult != 0) {
            // Jacobian计算失败，跳过此积分点
            continue;
        }

        // ====== 步骤5.2：计算Jacobian逆矩阵 ======
        double invJ[2][2];
        invJ[0][0] =  J[1][1] / detJ;
        invJ[0][1] = -J[0][1] / detJ;
        invJ[1][0] = -J[1][0] / detJ;
        invJ[1][1] =  J[0][0] / detJ;

        // ====== 步骤5.3：计算形函数导数（物理坐标系） ======
        double dN_dx[4], dN_dy[4];
        computeQuadShapeFunctionDerivatives(xi, eta, invJ, dN_dx, dN_dy);

        // ====== 步骤5.4：构造应变-位移矩阵B（3x8） ======
        double B[3][8];
        computeStrainDisplacementMatrix(dN_dx, dN_dy, B);

        // ====== 步骤5.5：获取材料切线刚度矩阵 ======
        double D[3][3];

        // 检查材料类型
        GPUMaterialType matType = static_cast<GPUMaterialType>(elem.materialType);

        if (matType == GPUMaterialType::PLASTIC_DAMAGE_CONCRETE_3D) {
            // 使用塑性损伤混凝土本构
            int stateIdx = elemIdx * 4 + gp_idx; // 假设每个单元4个积分点

            if (d_materialStates_committed != nullptr && d_materialStates_trial != nullptr) {
                // 准备材料参数
                PlasticDamageConcreteParams params;
                params.E = elem.materialProps[0];
                params.nu = elem.materialProps[1];
                params.ft = elem.materialProps[2];
                params.fc = elem.materialProps[3];
                params.beta = elem.materialProps[4];
                params.Ap = elem.materialProps[5];
                params.An = elem.materialProps[6];
                params.Bn = elem.materialProps[7];

                // 计算当前应变：ε_current = B * u
                double strainCurrent[6];
                double nu = elem.materialProps[1];  // 泊松比
                computeStrainFromBMatrix(B, elementDisplacements, strainCurrent, nu);

                // ====== 调试输出应变值 ======
                if (elemIdx == 0 && gp_idx == 0) {
                    printf("[GPU_DEBUG] elem=%d gp=%d: strain=[%.6e, %.6e, %.6e, %.6e, %.6e, %.6e]\n",
                           elemIdx, gp_idx,
                           strainCurrent[0], strainCurrent[1], strainCurrent[2],
                           strainCurrent[3], strainCurrent[4], strainCurrent[5]);
                    printf("[GPU_DEBUG] elem=%d gp=%d: nu=%.6e, E=%.6e\n",
                           elemIdx, gp_idx, nu, elem.materialProps[0]);
                }

                // 保存当前试算应变（用于commitState）
                if (d_strainCurrent_trial != nullptr) {
                    for (int i = 0; i < 6; i++) {
                        d_strainCurrent_trial[stateIdx * 6 + i] = strainCurrent[i];
                    }
                }

                // 读取历史应变（上次收敛的应变）
                double strainTotal[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
                if (d_strainTotal != nullptr) {
                    for (int i = 0; i < 6; i++) {
                        strainTotal[i] = d_strainTotal[stateIdx * 6 + i];
                    }
                }

                // 计算应变增量：Δε = ε_current - ε_history
                double strainIncr[6];
                for (int i = 0; i < 6; i++) {
                    strainIncr[i] = strainCurrent[i] - strainTotal[i];
                }

                // 调用混凝土本构（使用committed状态作为输入，计算trial状态）
                double stress[6];
                double D_tangent[36]; // 6x6刚度矩阵

                // ✅ 每次formTangent都从committed状态重新开始计算trial状态
                // 这与CPU材料的setTrialStrain()行为一致（见PlasticDamageConcrete3d.cpp:263-269）
                // CPU代码：eps_p = eps_pCommit; dp = dpCommit; dn = dnCommit; ...
                d_materialStates_trial[stateIdx] = d_materialStates_committed[stateIdx];

                computePlasticDamageConcrete3D(
                    strainIncr,
                    strainCurrent,    // 传递当前总应变而非历史应变
                    params,
                    d_materialStates_trial[stateIdx],  // 只修改trial状态，不碰committed
                    stress,
                    D_tangent
                );

                // 注意：不在formTangent期间更新应变历史！
                // strainTotal只在commitState()时更新，保存上次收敛的应变
                // 这样每次Newton迭代计算的strainIncr都是相对于上次收敛步的增量

                // ====== 调试输出材料状态和3D切线刚度 ======
                if (elemIdx == 0 && gp_idx == 0) {
                    printf("[GPU_DEBUG] elem=%d gp=%d: dp=%.4f dn=%.4f rp=%.3e rn=%.3e\n",
                           elemIdx, gp_idx,
                           d_materialStates_trial[stateIdx].dp,
                           d_materialStates_trial[stateIdx].dn,
                           d_materialStates_trial[stateIdx].rp,
                           d_materialStates_trial[stateIdx].rn);
                    printf("[GPU_DEBUG] elem=%d gp=%d: D_tangent[0,0]=%.3e D_tangent[2,2]=%.3e D_tangent[3,3]=%.3e\n",
                           elemIdx, gp_idx, D_tangent[0], D_tangent[14], D_tangent[21]);
                }

                // ✅ 使用静力凝聚从3D切线刚度(6x6)提取平面应力刚度(3x3)
                staticCondensation3Dto2D(D_tangent, D);

                // ====== 调试输出平面应力切线刚度矩阵 ======
                if (elemIdx == 0 && gp_idx == 0) {
                    printf("[GPU_DEBUG] elem=%d gp=%d: D_plane=\n", elemIdx, gp_idx);
                    printf("  [%.6e, %.6e, %.6e]\n", D[0][0], D[0][1], D[0][2]);
                    printf("  [%.6e, %.6e, %.6e]\n", D[1][0], D[1][1], D[1][2]);
                    printf("  [%.6e, %.6e, %.6e]\n", D[2][0], D[2][1], D[2][2]);
                }
            } else {
                // 如果没有状态数组，使用弹性刚度
                computePlaneStressD(E, nu, D);
            }
        } else {
            // 默认使用弹性材料矩阵
            computePlaneStressD(E, nu, D);
        }

        // ====== 步骤5.6：计算单元刚度矩阵增量 ======
        // K += B^T * D * B * detJ * t * w
        double factor = detJ * t * w;

        // ====== 调试：检查factor ======
        if (elemIdx == 0 && gp_idx == 0 && threadIdx.x == 0) {
            KERNEL_DEBUG_PRINT("[GPU_FACTOR] elem=%d gp=%d: detJ=%.6e t=%.6e w=%.6e factor=%.6e\n",
                   elemIdx, gp_idx, detJ, t, w, factor);
        }

        // 优化：先计算D*B（3x8矩阵）
        double DB[3][8];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 8; j++) {
                DB[i][j] = 0.0;
                for (int k = 0; k < 3; k++) {
                    DB[i][j] += D[i][k] * B[k][j];
                }
            }
        }

        // 然后计算B^T * (D*B)（8x8矩阵）
        for (int i = 0; i < 8; i++) {
            for (int j = 0; j < 8; j++) {
                double sum = 0.0;
                for (int k = 0; k < 3; k++) {
                    sum += B[k][i] * DB[k][j];
                }
                K_elem[i][j] += sum * factor;
            }
        }
    }

    // ====== 步骤6：将单元刚度矩阵写入全局内存 ======
    // ====== 调试输出单元刚度矩阵 ======
    if (elemIdx == 0) {
        printf("[GPU_DEBUG] elem=%d: K_elem diagonal = [%.6e, %.6e, %.6e, %.6e]\n",
               elemIdx, K_elem[0][0], K_elem[1][1], K_elem[2][2], K_elem[3][3]);
    }
    // 使用maxDOFPerElement确保与内存布局一致
    int offset = elemIdx * maxDOFPerElement * maxDOFPerElement;
    for (int i = 0; i < 8; i++) {
        for (int j = 0; j < 8; j++) {
            d_stiffness[offset + i * maxDOFPerElement + j] = K_elem[i][j];
        }
    }
}

/**
 * @brief GPU并行计算Truss单元刚度矩阵（简化实现）
 *
 * Truss单元：2节点，轴向受力
 * 刚度矩阵大小：6x6（2节点×3DOF）
 */
__global__ void computeTrussStiffness(
    GPUElementData* d_elementData,
    double* d_stiffness,
    double* d_nodeDisplacements,
    int numElements
)
{
    int elemIdx = blockIdx.x * blockDim.x + threadIdx.x;
    if (elemIdx >= numElements) return;

    GPUElementData& elem = d_elementData[elemIdx];

    // 检查单元类型
    GPUSupportedElementType elemType = getGPUElementType(elem.elementType);
    if (elemType != GPUSupportedElementType::TRUSS &&
        elemType != GPUSupportedElementType::TRUSS_SECTION &&
        elemType != GPUSupportedElementType::TRUSS2) {
        return;
    }

    // 提取节点坐标
    double x1 = elem.nodeCoords[0];
    double y1 = elem.nodeCoords[1];
    double z1 = elem.nodeCoords[2];
    double x2 = elem.nodeCoords[3];
    double y2 = elem.nodeCoords[4];
    double z2 = elem.nodeCoords[5];

    // 计算单元长度
    double dx = x2 - x1;
    double dy = y2 - y1;
    double dz = z2 - z1;
    double L = sqrt(dx*dx + dy*dy + dz*dz);

    if (L < 1e-14) return;  // 退化单元

    // 方向余弦
    double cx = dx / L;
    double cy = dy / L;
    double cz = dz / L;

    // 材料参数
    double E = elem.materialProps[0];  // 弹性模量
    double A = elem.materialProps[1];  // 截面积

    double k = E * A / L;

    // 局部刚度矩阵（6x6）
    double K_local[6][6];

    K_local[0][0] =  k*cx*cx; K_local[0][1] =  k*cx*cy; K_local[0][2] =  k*cx*cz;
    K_local[0][3] = -k*cx*cx; K_local[0][4] = -k*cx*cy; K_local[0][5] = -k*cx*cz;

    K_local[1][0] =  k*cy*cx; K_local[1][1] =  k*cy*cy; K_local[1][2] =  k*cy*cz;
    K_local[1][3] = -k*cy*cx; K_local[1][4] = -k*cy*cy; K_local[1][5] = -k*cy*cz;

    K_local[2][0] =  k*cz*cx; K_local[2][1] =  k*cz*cy; K_local[2][2] =  k*cz*cz;
    K_local[2][3] = -k*cz*cx; K_local[2][4] = -k*cz*cy; K_local[2][5] = -k*cz*cz;

    K_local[3][0] = -k*cx*cx; K_local[3][1] = -k*cx*cy; K_local[3][2] = -k*cx*cz;
    K_local[3][3] =  k*cx*cx; K_local[3][4] =  k*cx*cy; K_local[3][5] =  k*cx*cz;

    K_local[4][0] = -k*cy*cx; K_local[4][1] = -k*cy*cy; K_local[4][2] = -k*cy*cz;
    K_local[4][3] =  k*cy*cx; K_local[4][4] =  k*cy*cy; K_local[4][5] =  k*cy*cz;

    K_local[5][0] = -k*cz*cx; K_local[5][1] = -k*cz*cy; K_local[5][2] = -k*cz*cz;
    K_local[5][3] =  k*cz*cx; K_local[5][4] =  k*cz*cy; K_local[5][5] =  k*cz*cz;

    // 写入全局内存
    int offset = elemIdx * GPUElementConstants::TRUSS_STIFFNESS_SIZE;
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < 6; j++) {
            d_stiffness[offset + i*6 + j] = K_local[i][j];
        }
    }
}

/**
 * @brief GPU并行计算Beam单元刚度矩阵（TODO：待实现）
 */
__global__ void computeBeamStiffness(
    GPUElementData* d_elementData,
    double* d_stiffness,
    double* d_nodeDisplacements,
    int numElements
)
{
    // TODO: 实现Beam单元刚度矩阵计算
    // 包括2D梁和3D梁
}

//==============================================================================
// Phase 2 Optimization: GPU Node Coordinate Update Kernel
//==============================================================================

/**
 * @brief GPU并行更新节点坐标（Phase 2优化）
 *
 * 计算公式：currentCoord = initialCoord + displacement
 *
 * 优化策略：
 * - 避免CPU端遍历所有节点提取坐标
 * - 在GPU上直接计算当前坐标
 * - 预期收益：消除每步2-3ms CPU开销
 *
 * @param d_currentCoords 输出：当前节点坐标（GPU）
 * @param d_initialCoords 输入：初始节点坐标（GPU，首次缓存）
 * @param d_displacements 输入：节点位移（GPU，每步从CPU传输）
 * @param numNodes 节点数量
 */
__global__ void updateNodeCoordinates(
    double* d_currentCoords,        // 输出：当前坐标[numNodes*3]
    const double* d_initialCoords,  // 输入：初始坐标[numNodes*3]
    const double* d_displacements,  // 输入：节点位移[numNodes*3]
    int numNodes                    // 节点数量
)
{
    // 线程索引 = 节点索引
    int nodeIdx = blockIdx.x * blockDim.x + threadIdx.x;
    if (nodeIdx >= numNodes) return;

    // 更新3个自由度的坐标（x, y, z）
    for (int dim = 0; dim < 3; dim++) {
        int idx = nodeIdx * 3 + dim;
        d_currentCoords[idx] = d_initialCoords[idx] + d_displacements[idx];
    }
}

//==============================================================================
// End of Phase 2 Optimization
//==============================================================================

//==============================================================================
// C接口实现（用于从C++代码调用）
//==============================================================================

extern "C" {

/**
 * @brief 启动Quad单元刚度矩阵计算核函数
 */
int launchQuadStiffnessKernel(
    int numElements,
    const GPUElementData* d_elementData,
    double* d_stiffness,
    const double* d_nodeDisplacements,
    ConcreteState* d_materialStates_committed,
    ConcreteState* d_materialStates_trial,
    double* d_strainTotal,
    double* d_strainCurrent_trial,
    int maxDOFPerElement,
    cudaStream_t stream
)
{
    if (numElements <= 0) return -1;

    // 配置kernel启动参数
    int blockSize = 256;
    int gridSize = (numElements + blockSize - 1) / blockSize;

    // 启动kernel
    computeQuadStiffness<<<gridSize, blockSize, 0, stream>>>(
        const_cast<GPUElementData*>(d_elementData),
        d_stiffness,
        const_cast<double*>(d_nodeDisplacements),
        d_materialStates_committed,
        d_materialStates_trial,
        d_strainTotal,
        d_strainCurrent_trial,
        numElements,
        maxDOFPerElement
    );

    // 检查启动错误
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] launchQuadStiffnessKernel failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    return 0;
}

/**
 * @brief 启动Truss单元刚度矩阵计算核函数
 */
int launchTrussStiffnessKernel(
    int numElements,
    const GPUElementData* d_elementData,
    double* d_stiffness,
    const double* d_nodeDisplacements,
    cudaStream_t stream
)
{
    if (numElements <= 0) return -1;

    int blockSize = 256;
    int gridSize = (numElements + blockSize - 1) / blockSize;

    computeTrussStiffness<<<gridSize, blockSize, 0, stream>>>(
        const_cast<GPUElementData*>(d_elementData),
        d_stiffness,
        const_cast<double*>(d_nodeDisplacements),
        numElements
    );

    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] launchTrussStiffnessKernel failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    return 0;
}

/**
 * @brief 启动Beam单元刚度矩阵计算核函数
 */
int launchBeamStiffnessKernel(
    int numElements,
    const GPUElementData* d_elementData,
    double* d_stiffness,
    const double* d_nodeDisplacements,
    cudaStream_t stream
)
{
    if (numElements <= 0) return -1;

    int blockSize = 256;
    int gridSize = (numElements + blockSize - 1) / blockSize;

    computeBeamStiffness<<<gridSize, blockSize, 0, stream>>>(
        const_cast<GPUElementData*>(d_elementData),
        d_stiffness,
        const_cast<double*>(d_nodeDisplacements),
        numElements
    );

    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] launchBeamStiffnessKernel failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    return 0;
}

/**
 * @brief 启动节点坐标更新核函数（Phase 2优化）
 *
 * @param numNodes 节点数量
 * @param d_currentCoords 输出：当前节点坐标（GPU）
 * @param d_initialCoords 输入：初始节点坐标（GPU）
 * @param d_displacements 输入：节点位移（GPU）
 * @param stream CUDA流
 * @return 0-成功，-1-失败
 */
int launchUpdateNodeCoordinatesKernel(
    int numNodes,
    double* d_currentCoords,
    const double* d_initialCoords,
    const double* d_displacements,
    cudaStream_t stream
)
{
    if (numNodes <= 0) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] Invalid numNodes: %d\n", numNodes);
        return -1;
    }

    // 配置kernel启动参数
    int blockSize = 256;
    int gridSize = (numNodes + blockSize - 1) / blockSize;

    // 启动kernel
    updateNodeCoordinates<<<gridSize, blockSize, 0, stream>>>(
        d_currentCoords,
        d_initialCoords,
        d_displacements,
        numNodes
    );

    // 检查启动错误
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        KERNEL_DEBUG_PRINT("[GPU ERROR] launchUpdateNodeCoordinatesKernel failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    return 0;
}

}  // extern "C"
